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Abstract. A 3D dynamical model with a quasi- homogeneous core and a disk 
component is used for the chaos control in the central parts of elliptic;al galaxy. 
Numerical experiments in the 2D system show a very complicated phase plane 
with a large chaotic sea, considerable sticky layers and a large number of islands, 
produced by secondary resonances. When the mass of the disk increases, the 
chaotic regions decrease gradually, and, finally, a new phase plane with only 
regular orbits appears. This evolution indicates that disks in elliptical galaxies 
can act as the chaos controllers. Starting from the results obtained in the 2D 
system, we locate the regions in the phase space of the 3D system, producing 
regular and chaotic orbits. For this we introduce and use a new dynamical 
parameter, the S{w) spectrum, which proves to be useful as a fast indicator 
and allows us to distinguish the regular motion from chaos in the 3D potentials. 
Other methods for detecting chaos are also discussed. 
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1. INTRODUCTION 

Today it is believed, that dwarf elliptical galaxies seem to prefer the densest 
regions of the Universe and are found abundantly in galaxy clusters and groups 
(see BingeUi et al. 1987; Ferguson & Sandage 1989; Conselice et al. 2001). This 
circumstance most likely has important consequences for their evolution. Moore et 
al. (1998) have shown how late-type disk galaxies, that orbit in a cluster, can loose 
angular momentum by interactions with massive galaxies and, to a lesser degree, 
by tidal forces induced by the cluster potential. iV-body simulations performed 
by Mayer et al. (2001) show that small disk galaxies, that are close companions 
to a massive galaxy, will be affected likewise. A small disk galaxy is destabilized 
and develops a bar that gradually slows down, by dynamical friction, transporting 
angular momentum to the halo and to stars at larger radii. Since the latter are 
being stripped, angular momentum is lost. Gas is funneled in towards the center 
by torques exerted by the bar, where it is converted into stars, thus forming a 
nucleus. The small companion is heated by the subsequent buckling of the bar 
(see Mcrrifield & Kuijken 1999 and references therein) and by bending modes of the 
disk, and is transformed from a rotationally-flattened object into an anisotropic, 
slowly rotating spheroidal galaxy. The effect on a dwarf galaxy depends on its 
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orbit through the cluster around the massive companion. For instance, retrograde 
interactions have a much less damaging effect than progradc ones and may even 
preserve some of the initial disk structure. Thus, these simulations allow for the 
existence of fast-rotating dwarfs and for elliptical galaxies that still contain a stellar 
disk. 

During the last decade, dwarf elliptical galaxies, such as IC 3328, IC 0783, 
IC 3349, NGC 4431 and IC 3468, with spiral and barred structure have been 
discovered in the Virgo cluster (see Jcrjcn et al. 2000; Dc Rijckc ct al. 2001; 
Barazza et al. 2002; Simien & Prugniel 2002). Furthermore, Ryden et al. (1999) 
also report dwarf ellipticals with disky isophotes in the Virgo cluster, while De 
Rijcke et al. (2003) give photometric evidence for the presence of stellar disks in 
two Fornax dwarf galaxies FCC 204 and FCC 288. Moreover, there arc indications 
of the presence of disks in the giant elliptical galaxies NGC 83 and NGC 2320 (see 
Young 2002, 2005). The semi-analytic simulations of Khochfar & Burkert (2005) 
suggest that a few tens of percent of all disky elliptical galaxies, could have grown 
their stellar disks out of cold gas accreted from the intergalactic medium. A better 
understanding of the molecular gas in early-type galaxies could provide concrete 
evidence either for or against this disk scenario. 

Taking this into account, we considered interesting to construct a 3D dynamical 
model to study properties of the motion in an elliptical galaxy with a disk. In order 
to describe the motion in the region near the center of such a system, we use the 
potential 

V{x, y,z) = ^ [{x' + by' + cz') - ex {y' + z')] . 

V a;2 Jry'^Jr{a + VW^) 

(1) 

Potential (1) consists of two parts. The first term stands for the quasi-homogeneous 
core of the galaxy, the second term describes the disk (see Miyamoto & Nagai 
1975). The harmonic term of the potential is valid only within certain distances 
from the center of the galaxy. In our case, integration of orbits was done up 
to distances R — x"^ + y'^ + z'^ < 1, hinting that up to these distances from 
the galactic center, the harmonic term of Eq. (1) is valid. We must point out, 
that the perturbation term {—ex(y' + z'')) of the harmonic potential is an odd 
function with respect to x. We use this kind of perturbation term due to the fact 
that it has a finite energy of escape. Moreover, the first term of potential (1) is 
a deformed galactic model, which is approximately axisymmetric near the core 
but is deformed in its outer parts. A two-dimensional potential with this kind 
of perturbation was used by Contopoulos et al. (1987) in order to study large- 
scale stochasticity in a Hamiltonian system of two degrees of freedom, which may 
represent the inner parts of a deformed galactic model. In the present paper, the 
first term of potential (1), is the expansion of the potential used by Contopoulos et 
al. (1987) in a Hamiltonian system of three degrees of freedom. Similar potentials 
have been used in several papers, in order to study axisymmetric galactic models 
(see Contopoidos & Polymilis 1993; Siopis et al. 1995). 

In Eq. (1), b, c and e are parameters, is the mass, a and h are the scale 
length and the scale height of the disk, while oj is used for the consistency of 
the galactic units. We use a system of galactic units, where the unit of mass is 
2.325 X IO^Mq, the unit of length is 1 kpc and the unit of time is 0.997748 x 10® 
yr. The velocity unit is 10 km/s, while G is equal to unity. In the above units we 
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use the values: a; = 10 kms ^ kpc ^, b = 1.1, c = 1.5, e =1.08, a = 3, h = 0.125, 
while Md is treated as a parameter. 

The aim of this research is to investigate the role, played by the disk, on 
the regular or chaotic nature of orbits. Our investigation will be focused in the 
following: (i) we will try to connect the extent of the chaotic regions with the 
mass of the disk; (ii) we will look for sticky regions and secondary resonances 
introduced by the presence of the disk, and (iii) we will seek if there exist only 
one unified chaotic region, or different chaotic components. For this purpose, we 
shall use, apart of the classical methods, such as the Poincare phase plane and 
the Lyapunov Characteristic Exponents (LCEs) (Lichtenberg & Lieberman 1992), 
some modern methods such as the S{c) spectrum (Caranicolas & Papadopoulos 
2007; Caranicolas & Zotos 2010), the P{f) indicator (Karanis & Vozikis 2008) and 
a new dynamical parameter, the S{w) spectrum. 

The results are based on the numerical integration of the equations of motion 



X = 



y = -- 



d V{x,y,z) 

dx 

d V{x,y,z) 



dy 



d V{x, y, z) 

' = dz • 

The Hamiltonian to the potential (2) reads 

H=^{pl+pl+pl)+V{x,y,z)=E , (3) 

where Px,Py,Pz are the momenta per unit mass conjugate to x, y and z, while E 
is the numerical value of the Hamiltonian. 

The orbit calculations are based on the numerical integration of the equation 
of motion (2), which was made using a Bulirsh-Stocr routine in Fortran 95, with 
double precision in all subroutines. The accuracy of the calculations was checked 
by the consistency of the energy integral (3), which was conserved up to the twelfth 
significant figure. 

The present paper is organized as follows: in Section 2 wc present an analysis 
of the structure of the x — Px,y = 0,Py > Poincare phase plane of the 2D system 
and the different families of orbits. In the same section a study of the evolution of 
the sticky regions is presented. Moreover, we investigate the evolution of different 
chaotic components of the system. In Section 3 we study the character of orbits 
in the 3D system, using a new dynamical indicator, the S(w) spectrum. Special 
interest is given to the evolution of the sticky regions and the chaotic components. 
In Section 4, a discussion and the conclusions of this research are presented. 



2. THE STRUCTURE OF THE x - PHASE PLANE: CHAOTIC 
COMPONENTS AND STICKY REGIONS 

In this section we will analyze the structure of the a;— 5 2/ = OjPj/ > Poincare 
phase plane of the corresponding 2D potential, i.e., we will consider orbits on the 
galactic plane (2: = 0). The corresponding Hamiltonian is 



H2 = l{pl+pl)+V{x,y)=E2 , 



(4) 
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where E2 is the numerical value of the Hamiltonian. 

Figures 1 (a-d) show the x — Px phase plane of the 2D system. Let us start 
from Figure la, where Md = 100, E2 = 20. One observes a phase plane with a 
complicated structure. Several chaotic components seem to exist. There are also 
two main regions of regular motion. The first regular region, forms invariant curves 
around the two invariant points on the x-axis, while the second region is composed 
of two sets of islands intersecting the ^^.-axis. Both the above islands are produced 
by quasi-periodic orbits characteristic of the 1:1 resonance. A considerable regular 
region is occupied by two additional sets of three islands symmetrical with respect 
to the X-axis. Furthermore, there are sets of smaller islands embedded in the 
chaotic sea and produced by secondary resonances. Several sticky regions seem to 
be present, but we shall come to this interesting point later in this Section. Figure 
lb shows the phase plane when = 200, E2 = —10. At a first glance, we see 
that the regular region has increased, while the chaotic regions have decreased. 
Note that the two invariant points near the px-Bxis are now stable. The chaotic 
components are present in this case too. The different chaotic components will be 
studied later in this Section. We must also point out, that in this case we can see 
one main sticky region and some smaller secondary sticky regions as well. Figure 
Ic shows the x — Px phase plane in the case when M4 = 600, E2 = —132. Here 
things are different. Almost the entire phase plane is regular, while only a small 
chaotic layer with some tiny sticky regions is observed. On the other hand, all 
periodic points are now stable, and no secondary resonances are observed. Figure 
Id is similar to Figure Ic, but for Md — 1200, E2 = —315. In this case, the entire 
phase plane is covered by regular orbits, while the percentage of chaotic orbits, if 




Fig. 2. Orbits in the 2D potential. The initial conditions in panel (a): xo — —0.90, 
Pxo = 0, panel (b): xo — 0.32, pxo ~ 5.2, panel (c); xo — 0.05, PxO = 6.2, panel (d): 
xo = —0.952, pxo = 0, panel (e): xo = —0.01, Pxo = 4.7, panel (f): xo = —0.01, PxO = 6, 
panel (g): xo — 0.327, Pxo — 1.16, panel (h): a::o ~ 0.09, Pxo = 0. In panels (a-e) 
Md = 100, E2 = 20, in panels (f-h) Ma = 200, E2 = -10. 



any, is negligible. 

Figures 2 (a li) show several representative orbits in the 2D system. In all 
orbits yo = 0, while the value oipyo is always found from the energy integral. The 
values of initial conditions and parameters are given in the caption. All orbits are 
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Fig. 3. A plot of the chaotic percentage A% vs. Ma- 



regular except of the orbit given in Figure 2h, which is chaotic. The integration 
time for all 2D orbits shown in Figures 2 (a-h), is 100 time units. 

Thus we conclude, that in elliptical galaxies with massive disks in the central 
regions a decrease of chaos is expected, while the percentage of chaotic orbits is 
larger, when a small disk is present. But the most important conclusion drawn 
from the above numerical study, is that disks in the centers of elliptical galaxies 
can act as chaos controllers. Figure 3 shows the percentage A% of the surface of 
section occupied by chaotic orbits vs. Md- We see that A% tends asymptotically 
to zero, when the mass of the disk increases. 

Now let us study the evolution of the sticky orbits of the 2D system. We 
shall follow in detail the evolution of a sticky orbit with the initial conditions: 
Xq = —0.962, Pjio = 0, the values of all other parameters arc as in Figure la. The 
results are given in Figure 4 (a-d). Figure 4a shows the sticky region formed in 
the X — Px phase plane for a time period of about 1200 time units. After that 
the test particle leaves the above sticky region, entering to a larger sticky region 
shown in Figure 4b. There it stays to about 75 000 time units and then it moves 
to the third sticky region shown in Figure 4c, where it stays until 3 x lO*^ time 
units. After reaching this stage, our numerical calculations were not continued. 
Our feeling is that the evolution of the sticky orbit was completed here. Actually 
we observe a hierarchy regarding the sticky regions. We believe that in Figure 4c 
we see a chaotic component of the 2D system, not a sticky region. Figure 4d shows 
a plot of the LCE of the sticky orbit for a time period of 3 x 10^ time units. 

A better view of the sticky orbit evolution can be seen using the S{c) spectrum. 
Figure 5a shows the 5(c) spectrum for an orbit producing a set of eight small 
islands, shown in Figure la. The initial conditions are: xq = —0.952, yo = Pxa = 0, 
and the values for all other parameters are as in Figure la. As expected, we 
observe eight well defined [/-type spectra. The motion is regular. Figure 5b shows 
the S{c) spectrum for the sticky orbit. Here, we can see again eight spectra, each 
corresponding to an island. The basic difference between Figures 5a and 5b is that 
in Figure 5b wc observe a large number of asymmetric peaks. Those additional 
peaks indicate the sticky motion. It is well known, that in the dynamical systems 
of two degrees of freedom, sticky orbits are the orbits which stay for long time 
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Fig. 4. Panels (a-c): evolution of a sticky orbit; panel (d): LCE vs. time for the 
sticky orbit. 

periods near the last Kolmogorov-Arnold-Moser (hereafter KAM) torus, before 
they escape to the surrounding chaotic sea (see Karanis & Caranicolas 2002). 

According to the Kolmogorov-Arnold-Moser theorem, most orbits lie on tori in 
the phase space. However, near all unstable periodic points, there is some degree 
of stochasticity. This is best seen on the surface of section, on which the successive 
intersections of orbits passing close to the unstable periodic points, do not lie, in 
general, on closed invariant curves, but fill stochastically a certain defined area. 
These orbits are called stochastic, or chaotic, or semiergodic. The phenomenon 
of the onset of large scale chaoticity, has been studied in more detail in recent 
years. It was found, that the invariant curves that separate chaotic regions in 
the neighborhood of two unstable periodic orbits are destroyed as the energy goes 
beyond a critical value. Thus, we have communication between the two chaotic 
regions. The critical value of the energy occurs when the last KAM curve (or 
the last KAM torus in the phase plane) separating the two resonant regions is 
destroyed. 

Now let us go to the evolution of the sticky orbit, using the S{c) spectrum. 
During the first sticky period, which is about 1100 time units, in Figure 5b one 
observes eight separate complicated spectra. In Figure 5c the time is 1500 time 
units, and the eight spectra are very similar to those seen in Figure 5b. Note 
that here the eight spectra are connected. This indicates that the test particle 
has left the first sticky region in order to continue its wandering in a larger sticky 
region. Figure 5d, shows the S{c) spectrum for 70 000 time units. We see that the 
spectrum tends to take characteristics of a chaotic spectrum. Note that the time 
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Fig. 5. Panel (a): the S{c) spectrum for an orbit producing tlie set of eight small 
islands, shown in Fig. la; panels (b-d): evolution of the S{c) spectrum of a sticky orbit. 
Details are given in the text. 

intervals of the sticky periods, obtained by the evolution of the S{c) spectrum, are 
very close to those obtained by the formation of the x — Px phase planes, shown 
in Figxnes 4 (a d). 

The next step is to investigate different chaotic components of the 2D system. 
The corresponding results are shown in Figures 6 (a-d). Figure 6a shows the 
chaotic region formed by an orbit with the initial conditions: xq = 0.0, Pxo = 4.6 
over a time period of 3 x 10'^ time units. This is the second chaotic component 
of the 2D system. Figure 6b shows the third chaotic component which is formed 
by an orbit with the initial conditions: xq = —0.51, PxO = 6.0. The time period is 
3 X 10^ time units. Figure 6c shows all the chaotic components together. The LCEs 
of the above three chaotic components arc shown in Figure 6d. The numbers 1,2 
and 3 correspond to the first, second and third chaotic component, respectively. 
As we see, each chaotic component has a different value of LCE (see Saito & 
Ichimura 1979). Here we must note that a hierarchical structure in the 2D model 
is displayed not only by the stickiness but also by chaos. 

3. ORDER AND CHAOS IN THE 3D SYSTEM 

Now let we switch to study the behavior of the 3D system. We take the initial 
conditions {xo,Pxo, zo),yo = PzO = 0, where {xo,Pxo) is a point on the phase planes 
of the 2D system. Obviously this point lies inside the limiting curve 

^-pl + Vix)=E, , (5) 

which is the curve containing all the invariant curves of the 2D system. We take 
E = E2, and the value of Pyo for all orbits is obtained from the energy integral 
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Fig. 6. Panel (a): chaotic component 2; panel (b): chaotic component 3; panel (c;. an 
the chaotic components together; panel (d): the LCEs for the three chaotic components 
of the 2D system. 



(3). Our numerical calculations in the 3D system show, that the orbits with the 
initial conditions {xq, p^o, zq), uq = pzo = 0, where (xq, Pxo) is a point in the 
chaotic regions of Figures 1 (a-d), for all permissible values of Zq give chaotic 
orbits. On the other hand, it would be interesting to know what happens with 
different chaotic components observed in the 2D system, shown in Figures 6 (a-d). 
The question is, do they merge to form a unified chaotic region in the 3D space, 
or they continue to exist as three separate chaotic components in the 3D system? 
In order to give an answer, we computed the LCEs for three orbits, each starting 
in one of the three different chaotic components, for a time period of 10^ time 
units. The results are shown in Figure 7. The numbers 1, 2 and 3 indicate the 
three chaotic components which have different LCEs values. This result seems to 
be consistent with the outcomes obtained by Cincotta et al. (2006), where in a 
3D system with divided phase space, separate chaotic components actually exist. 

One may ask an interesting question: what is the nature of orbits which have 
initial conditions (a;o, Pxo)^ zq, yo = Pzo = 0, where (a;o, Pxo) is a point in each 
regular region of Figures 1 (a-d)? In order to give an answer, we will introduce 
and use a new type of dynamical indicator, the S'(w) spectrum. The parameter 
Wi is defined as 



{Xj -Pxi) - jZj ~Pzi) 
Pyt 



(6) 



where (xi, Zi, Pxi, Pyi, Pzi) are the successive values of {x, z,Px,Py,Pz) of the 3D 




Fig. 8. Panel (a): a 3D regular orbit; panel (b): a plot of the LCE vs. time for the 
orbit shown in (a); panel (c): the 3(111) spectrum of the orbit shown in (a); panel (d): 
the P{f) indicator for the orbit shown in (a). See the text for details. 
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Fig. 9. The same as in Figures 8 (a-d) but for a resonant 3D orbit. See the text for 
details. 



orbit. The dynamical spectrum of the parameter w is its distribution function 

AN{w) 



Siw) = 



NAw 



(7) 



where AN{w) is the number of the parameters w in the interval w,w + Aw after 
N iterations. In order to study the character of a 3D orbit, the 5(c) spectrum can 
be also used. Note that the coupling of the third component, z, carrying all the 
information about the 3D motion, is hidden in the definition of the S{c) spectrum, 
but in any case it affects the values of x,pj. and py. Using the definition of the 
S{w) spectrum, we overtake this minor drawback and create a new dynamical 
spectrum, suitable especially for 3D orbits. 

Figures 8 (a d) show the results for a 3D regular orbit. The orbit which is 
shown in Figure 8a has the initial conditions: xq = 0.82, yo = p^o = PzO = 0, 
Zq = 0.01, while for all 3D orbits the value oipyQ is always found from the energy 
integral (3). The corresponding values of all other parameters are as in Figure la. 
The LCE of the orbit, shown in Figure 8b, vanishes indicating the regular motion. 
Figure 8c shows the S{w) spectrum of the orbit. This is a well defined [/-type 
spectrum characteristic for the regular motion. In Figure 8d, we see the P{f) 
indicator which also indicates the regular motion. Figures 9 (a-d) are similar to 
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Fig. 10. The same as in Figures 9 (a-d) but for a resonant 3D orbit with different 
initial conditions. See the text for details. 



Figures 8 (a-d) but correspond to a 3D resonant orbit with the initial conditions: 
Xq = 0.27, 2/0 = Pzo = PxQ = 5-3, zq = 0.01, producing three f7-type spectra. 
The values of all other parameters arc as in Figure lb. 

Figures 10 (a-d) are similar to Figures 9 (a-d) but for a 3D resonant orbit 
with the initial conditions: Xq = 0.47, j/o = PzO = 0, PxO = 4.9, Zq = 0.01, 
producing five ?7-type spectra. The values of all other parameters are as in Figure 
lb. Finally, in Figure 11 (a d), wc present the results for a 3D chaotic orbit. The 
initial conditions are: xq = 0.09, j/o = Pxa — PzO = 0) -^o = 0.1. The values of 
all other parameters are as in Figure la. Once more, we observe that the results 
from all dynamical indicators, regarding the character of motion, coincide. The 
integration time for all 3D orbits is 100 time units, for the S{w) spectrum it is 
2 X 10'* time units and for the P{f) indicator it is 4 x 10^ time units. 

Using the above method, we have computed a large number of orbits (about 
1000) in the 3D dynamical system. Numerical outcomes suggest that for all {xq, 
Pxo) in the regular regions of Figures 1 (a d) and for small values of zq (zq < 0.12) 
the motion is regular, while for all larger permissible values of zq the motion 
becomes chaotic. 
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Fig. 11. The same as in Figures 8 (a-d) but for a chaotic 3D orbit. See the text for 
details. 



4. DISCUSSION AND CONCLUSIONS 

In this article we have studied the properties of motion in a 3D Hamiltonian 

system, describing the motion in the inner parts of a deformed galactic model. 
We started our investigation from the 2D system, because orbits confined in the 
galactic plane {z = 0) display some very interesting features, such as sticky regions, 
chaotic components and islandic motion produced by a large number of secondary 
resonances. 

Results from the study of the phase planes indicate, that there is an hierarchy 
in sticky regions. A test particle can stay in sticky regions for a time period 
of about 75 000 time units, before leaving to the corresponding chaotic region. 
Several chaotic components are also observed in the 2D system, each one having 
a different value of LCE. Thus we conclude, that both the stickiness and chaos 
display a hierarchical structure in the 2D model. An interesting result of this 
investigation is that the percentage A% of the surface of the section occupied by 
chaotic orbits decreases, tending asymptotically to zero, when the mass of the 
disk increases. This suggests that disks in elliptical galaxies can act as the chaos 
controllers. 

Our numerical calculations in the 3D system show, that orbits with the initial 
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conditions (xq, PxOi zq), yo = Pzo = 0, where {xq, Pxo) is a point in the chaotic 
regions of Figure 1 (a-d), for all permissible values of zq give chaotic orbits. Using 
the new S{w) spectrum we find that orbits with the initial conditions {xq, PxO^ 
Zq), yo = PzO = 0; where (a;o, Pxo) is a point in the regular regions of Figure 1 
(a-d), for small values of zq are regular, while for larger values of zq they become 
chaotic. It is also of particular interest that the 3D system displays three different 
chaotic components, not a unified chaotic region. 

A very effective and reliable tool to distinguish between the regular and chaotic 
motion in the 3D dynamical systems is the new S{'w) spectrum. This spectrum, 
which is an advanced form of the 5(c) spectrum, allows us to detect islandic 
3D motion of the resonant orbits, since it produces as much spectra as the total 
number of islands in the x — px, y — Q, Py > surface of section (see Figures 
9a and 10a). Moreover, the comparison with other dynamical parameters, such 
as the LCE and the P{f) indicator, shows that the results obtained by the S{w) 
spectrum are reliable. 
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